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Negative static permittivity and violation of Kramers-Kronig relations in 

qnasi-two-dimensional crystals 
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We investigate the wave-vector and frequency-dependent screening of the electric field in atomi¬ 
cally thin (quasi-two-dimensional) crystals. For graphene and hexagonal boron nitride we find that, 
above a critical wave-vector q^, the static permittivity e{q > qc,ijj = 0) becomes negative and the 
Kramers-Kronig relations do not hold for e{q > qc,a}). Thus, in quasi-two-dimensional crystals, 
we reveal the physical confirmation of a proposition put forward decades ago (Kirzhnits, 1976), 
allowing for the breakdown of Kramers-Kronig relations and for the negative static permittivity. 

In the vicinity of the critical wave-vector, we find a giant growth of the permittivity. Our results, 
obtained in the ab initio calculations using both the random-phase approximation and the adiabatic 
time-dependent local-density approximation, and further confirmed with a simple slab model, allow 
us to argue that the above properties, being exceptional in the three-dimensional case, are common 
to quasi-two-dimensional systems. 

PACS numbers: 77.22.Ch, 73.22.Pr 


The concept of causality plays one of the central roles 
in contemporary science [1]. It is well known that causal¬ 
ity in the time-domain (the impossibility for an effect 
to precede the cause in time) leads to the analyticity 
of a causal response-function in a complex half-plane in 
the frequency-domain, which, in turn, leads to Kramers- 
Kronig (KK) relations between the real and the imagi¬ 
nary parts of the response function [2]. 

It must be, however, recognized that the causality as¬ 
sumes that the response-function is applied to a cause 
and it produces an effect. In the case of the longitudinal 
electric field in a translationally invarient or a periodic 
system, the definition of the permittivity £(q, w) reads 
0tot(q,w) = (()ext(q,w)/e(q,a;), where (()ext and (ptot are 
the scalar potentials of the externally applied and the 
total electric fields, respectively. Since the cause is ^ext 
and the effect is (ftot, not vice versa, this is 1/e that is 
guaranteed to be causal, but not e itself [3]. Accordingly, 
KK relations must be satisfied by 1/e, but may or may 
not be satisfied by e. For |q| > 0, this leaves e(q, w = 0) 
a freedom to be negative without violating the causality 
or destroying the stability of the system [4, 5]. If this 
happens, then the inverse permittivity has zeros in the 
upper half of the complex w-plane, making the permit¬ 
tivity itself a non-analytic function. 
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where xgg' {z, z', q, w) is the density-response function of 
the system in the mixed, reciprocal in the system plane 
(xy) and real in the z-direction representation (G are the 
2D reciprocal lattice vectors). 
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FIG. 1. (Color online) Schematics of 2D material under an 
external field, a) Q2D single-layer geometry and b) 3D super¬ 
cell geometry. 


In the three-dimensional world the realizations of such 
negative static permittivity are scarce and they mostly 
concern exotic non-crystalline systems [6-10]. In this 
work we show that, above a critical wave-vector q > 
in the first Brillouin zone, the permittivity e(q, w) of the 
quasi-two-dimensional (Q2D) systems of the monolayer 
graphene and boron nitride is negative in the static limit. 
Accordingly, KK relations for the permittivity do not 
hold in this case. The inverse permittivity, on the con¬ 
trary, remains causal and does satisfy KK relations. 

We start by writing the permittivity of a Q2D crystal 
[11] (atomic units e^ = h = rrif. = 1 are used throughout 


Our time-dependent density-functional theory 
(TDDFT) calculation of the permittivity consists of 
two steps. Since Q2D systems lack periodicity in the 
z-direction, it is customary to use the super-cell method 
[12-14]. First, in the super-cell geometry, we calculate 
the density-response function XGg.G'g'(Q j ^5 of 
auxiliary 3D system comprised of an infinite periodic 
array of monolayers with the separation d between them, 
as is schematized in Fig. 1 b) (the dependence of x on d 
is shown explicitly). Here g are reciprocal vectors in the 
z-direction. 

We have conducted the calculation for the monolayer 
pristine graphene using the full-potential linear aug- 
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merited plane-wave (FP-LAPW) code Elk [15]. The 2 - 
axis period d of the super-cell was taken 20 a.u. The 
fc-point grid of 512 x 512 x 1, 30 empty bands, and the 
damping parameter of 0.002 a.u. were used in both the 
ground-state and the linear-response calculations. The 
former was carried out within the local-density approxi¬ 
mation (LDA) [16] for the exchange-correlation (xc) po¬ 
tential, while the latter was the random-phase approxi¬ 
mation (RPA) one (i.e., the xc kernel fxc [17] was set to 
zero). 

Results for £ 3 ^, obtained through 
1 Att 

-7-n ^ —TXoo,oo(q, w; d), ( 2 ) 

£3D(q,w;d) 


are presented in the left panels of Figs. 2 and 3 for 
q = 0.049 and 0.152 a.u., respectively, along the TM 
direction. 

It is, however, known that e 3 D(q,w;d), calculated in 
the super-cell geometry, is a quantity completely differ¬ 
ent from the permittivity £(q, w) of a single layer [ 11 - 
14], as can be also immediately appreciated from the 
d-dependence of the former. Our second step consists, 
therefore, in finding the density-response function y of 
the single-layer system from that of the array of those 
layers y. This can be conveniently done by virtue of the 
matrix relation [ 11 ] 

X(q,w) = y(q,a;) [1 -h C'(q)y(q, w)]"^, (3) 


where the elements of the matrix C are given by 
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FIG. 2. (Color online) Left: 3D permittivity of the array 
of graphene layers. Right: Permittivity of a single graphene 
layer. The wave-vector q = 0.049 a.u. is below the critical 
value Qc ~ 0.118 a.u. 



In particular, y calculated by Eqs. (3)-(4) is free of the 
spurious inter-layer interaction, which is present in y. 

By the use of Eqs. (3)-(4), we find y in the 3D 
reciprocal-space representation. Then, by the inverse 
Eourier transform to the mixed representation and using 
Eq. (I), we obtain the permittivity £(q, w). The latter is 
plotted in the right panels of Figs. 2 and 3. 

A striking feature in Fig. 3, right panel, is that £(q, w = 
0 ) is negative and, since Ime(q, w > 0 ) > 0 , the per¬ 
mittivity of graphene does not satisfy Kramers-Kronig 
relations [ 2 ] 
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V denoting the principal value of the integrals, with 
R{u}) = £(q,a;). This fact is further illustrated in Fig. 4, 
where the real part of £(q, a;) is compared with the KK 


FIG. 3. (Color online) Left: 3D permittivity of the array 
of graphene layers. Right: Permittivity of a single graphene 
layer. The wave-vector q — 0.152 a.u. is above the critical 
value qc « 0.118 a.u. 


transform of its imaginary part: The two functions coin¬ 
cide atq<qc (left panel), but they are largely different at 
q > qc (right panel). We found the critical wave-vector 
for graphene to be qc ~ 0.118 a.u. (0.223 A“^). [18] 
On the other hand, it can be seen in Figs. 2 and 3, left 
panels, that the array system has a positive static per¬ 
mittivity, which cannot be otherwise for a 3D periodic 
system within RPA [5]. 

Since KK relations are not satisfied by £(q,a;), the lat¬ 
ter must have a singularity in the complex uj upper half¬ 
plane. The singularity can only be a pole at a; = a;s 
satisfying 


1 

e(q,Ws) 


= 0 . 


( 6 ) 


Considering that (a) l/£(q, a;) is a real continuous func- 
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FIG. 4. (Color online) Real and imaginary parts of the per¬ 
mittivity of a single-layer graphene. Separately is plotted the 
real part obtained from the imaginary part by the use of KK 
relation. Left: q < q^, KK relation holds. Right: q > q^, KK 
relation does not hold. 


tion on the positive imaginary axis of the w-plane; (b) 
l/e(q,a; = 0) < 0 for g > qd and (c) l/e(q, w = too) = 1, 
we conclude that, at q > qc, there exists a point lUs on 
the positive imaginary axis of the w-plane which satisfies 
Eq. (6). To find this point, we write by virtue of Cauchy’s 
integral formula 
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Expanding the complex inverse pemittivity in the right- 
hand side of Eq. (7) via its real and imaginary parts and 
using the parity properties of those functions, we can 
write 
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FIG. 5. (Color online) Inverse permittivity of a single 
graphene layer as a function of the imaginary frequency. At 
q < qc (black dashed curve), there is no zero (the permittiv¬ 
ity is analytic in the upper complex tj-plane). At q > qc (red 
solid curve), the inverse permittivity has a zero (indicated 
by a circle). Accordingly, the permittivity has a pole at this 
frequency. 


happen below the critical wave-vector. In Appendix A, 
we compare our results with the analytical ones known 
in the low-g regime [19, 20]. 

We can gain a further insight into the situation us¬ 
ing the approximate analytical relation between the Q2D 
permittivity of a single layer and the 3D permittivity of 
the array of those layers 
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rather than with the ’exact’ numerical solution of 
Eqs. (3)-(4). Equation (10), derived in Ref. 14, is a good 
approximation at q far from the critical value from the 
both sides, as we demonstrate below in Fig. 7. Solving 
Eq. (10) with respect to £ 30 , we find that 1/e is zero if 


£3D(q,a;;d) = 1-- 

qd 
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Further, the equality of the second term on the right- 
hand side to the third one can be easily proven with the 
use of the KK relations for l/£(q, w). We then have 
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We use Eq. (9) to calculate the inverse permittivity 
on the positive imaginary w-axis from our results for it 
on the real axis. In Fig. 5, this is plotted for the two 
wave-vectors, below and above the critical value. Above 
the critical wave-vector, the inverse permittivity crosses 
zero (indicated by a circle in Fig. 5), which does not 


In Fig. 6 , we plot £ 30(^1 w; d) along the positive imag¬ 
inary w. Although £ 3 D(q, w;d) is analytic in the upper 
complex w-plane at all values of q, it gives rise to a zero 
in l/£ (a pole in e) when the condition (11) is met. In 
Fig. 6 this is shown as an intersection, in the case of 
q > qc, with the straight horizontal line representing the 
right-hand side of Eq. (11). 

Figure 7 is presented in support of the fact that the 
permittivity obtained through the ’exact’ numerical pro¬ 
cedure via Eqs. (3), (4), and (1) can be accurately ap¬ 
proximated by the simple analytical formula of Eq. (10), 
if q is sufficiently below or above the critical wave-vector 
(upper panels). On the contrary, the same comparison 
done for the wave-vector slightly below and above the 
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FIG. 6 . (Color online) The permittivity of an array of 
graphene layers as a function of the imaginary frequency. 
Horizontal dashed lines show the values of £30 at which, by 
Eq. (11), the permittivity of a single-layer graphene may be¬ 
come singular. This never happens q < (black dashed 
curve), while this happens at q > q^ (red solid curve). 
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FIG. 7. (Color online) The permittivity of single-layer 
graphene calculated with the ’exact’ numerical procedure us¬ 
ing Eqs. (3), (4), and (1), and the approximate analytical 
Eq. (10). Upper panels: The wave-vector is well below (left) 
and above (right) the critical value qc ~ 0.118 a.u. Lower pan¬ 
els: The wave-vector is slightly below (left) and above (right) 

Qc- 


critical value (lower panels), reveals the complete inap¬ 
plicability of the approximate formula (10) in the vicinity 
of the critical wave-vector. Moreover, a giant increase in 
the absolute value of the permittivity occurs close to the 
critical wave-vector. 

Graphene is known to be a semi-metal, possessing a 
remarkable property of Dirac’s cones touching in the K- 
point of its band-structure [21]. A natural question arises 
whether the negative static permittivity and the violation 
of KK relations in graphene are in any way related to the 
Dirac’s cones in this material. To answer this, in Fig. 8 


FIG. 8 . (Golor online) The permittivity of 2D hexagonal 
boron nitride below (left) and above (right) the critical wave- 
vector qc ~ 0.323 a.u. 


we present results for the permittivity of the hexagonal 
boron nitride (h-BN), known to be an insulator [22]. Sim¬ 
ilar to graphene, above a critical wave-vector r; 0.323 
a.u. (0.610 A-i) (right panel of Fig. 8), the permittiv¬ 
ity of BN does not satisfy KK relations, while having 
a negative static limit. Furthermore, in Appendix B, 
we demonstrate that a simple local model of a metallic 
slab in vaccum supports the negative static permittiv¬ 
ity at larger wave-vectors. This shows that the negative 
static permittivity and the breakdown of KK relations is 
a rather general property common to Q2D systems. 

Importantly, in perfect 3D crystals the negative static 
permittivity is only possible due to the dynamic xc effects 
in the electronic response [5]. Since our results for Q2D 
crystals are obtained within the RPA, i.e., neglecting the 
xc effects, and the negative static permittivity occurs 
possible, the situation is fundamentally different with 
Q2D crystals: This is the finite but microscopic thickness 
of the crystal, which is also the break of the periodicity in 
one dimension, that makes the negative static permittiv¬ 
ity possible. Nonetheless, static and dynamic many-body 
effects play an important part in Q2D crystals [23, 24], 
which have not been accounted for in the present study. 
In Appendix C, we show that the the inclusion of the xc 
kernel fxc on the level of the adiabatic time-dependent 
local-density approximation (ATDLDA) [17, 25, 26] does 
not lead to a significant change in the results. The inclu¬ 
sion of the same effects within TDDFT with more elab¬ 
orate fxc, e.g., following the schemes known in the 3D 
case [27, 28], presents a challenge in the case of Q2D sys¬ 
tems. Furthermore, for Q2D crystals supported on sub¬ 
strates, the interaction with the latter strongly influences 
the excitation processes [29], which is also a demanding 
problem to be addressed in the future. 

For the accurate interpretation of the results, it is nec¬ 
essary to keep in mind the exact meaning of the permit¬ 
tivity (1) of a Q2D crystal. This definition is given in two 
steps: [11] First, the 2D conductivity (t^d respect to 
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the external field is introduced 

j2D(q,w) = 0-20 (12) 


where w) is the uniform in the ^-direction exter¬ 

nal electric field, and j 2 D(q, w) is the 3D current-density 
integrated in the z-direction and averaged over the unit 
cell in the xy-plane. Secondly, the permittivity of a Q2D 
crystal is defined by the relation 


£(q,w) 


= 1 + 


"2D 


(q,w). 


(13) 


rigorously valid for a strictly 2D system. The final justi¬ 
fication of Eq. (13) is that with this definition the usual 
formula for the energy dissipation 

(5(q,w) =-^|Eext(q,w)plni—(14) 
Airq e(q,w) 


holds for a Q2D crystal exactly. However, as detailed in 
Ref. 11, the Q2D permittivity cannot be attributed the 
meaning of the coefficient of proportionality between the 
external and the total fields and, hence, Eq. (14) cannot 
be rewritten in terms of Etot and Ime. We note that 
these complications call for the particular caution in the 
consideration of the negative static permittivity in the 
context of the 2D superconductivity [4, 5]. The behav¬ 
ior of the total field along the z-direction in graphene is 
further discussed in Appendix Sec. D. 


In conclusions, we have established the violation of 
Kramers-Kronig relations by the wave-vector and fre¬ 
quency dependent permittivity of quasi-two-dimensional 
graphene and boron nitride above a critical magnitude 
of the wave-vector, and the static permittivity was found 
negative in this case. The mechanism for the negative 
static permittivity was shown conceptually different from 
that in the 3D case: It is due to the system finite micro¬ 
scopic thickness rather than to the exchange-correlation 
effects. Our findings suggest the fundamental differences 
of the screening and the electronic excitation processes in 
quasi-2D crystals as compared with both 3D and purely 
2D systems. It is, however, discussed that further work 
is required to consider the present results in the context 
of the 2D superconductivity. 
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Appendix A: Low-Q permittivity of graphene 



u (eV) 


FIG. 9. Inverse permittivity of graphene as a function of the imaginary frequency calculated within the framework of this paper 
(solid lines) and its analytical long-wave behaviour (Al) within a 2D model [19, 20] (dashed lines). 


At small q, the permittivity of graphene, thought as a strictly 2D system, can be written as a function of the 
imaginary frequency as [19, 20] (the density-response function of Eq. (2) in Ref. 19 must be two times less) 


s{q, iu) = 1 + 


irq 

2 \/ m ^ + {qvf)'^' 


(Al) 


where Vf ^ 0.26 a.u. [19] is the Fermi velocity. In Fig. 9 we compare our first-principles results for graphene as a 
Q2D system to those from Eq. (Al). The conclusions are as follows: 

1. At very small wave-vector (q — 0.003 a.u.) the two calculations are in good agreement except for u at and very 
close to zero. The latter disagreement at very small q and u (or w, in the real-frequency calculation) is due to a 
the principle difficulty in achieving a good accuracy in the numerical calculation with finite fc-grid and damping 
in the very vicinity of the Dirac point. In this regime results obtained through Eq. (Al) can be thought superior 
to the ab initio ones. See also discussion of Fig. 6 in Ref. 11. 

2. Small but not too small wave-vector {q = 0.009 and 0.021 a.u.). The agreement between the two calculations is 
very good in the low u range. There is, however, no reason for them to agree at higher u, since the analytical 
formula is built on the two-bands model, while the numerical calculation uses the realistic band-structure (30 
bands for this calculation) (see also Fig. 6 in Ref. 11 to illustrate the same point in the real-frequency calculation). 

3. Larger wave-vector {q = 0.049 a.u.). The two calculations disagree, the small-q expansion being not relevant 
any more. 

In both regimes 2 and 3 the numerical results are superior to the analytical ones. 


Appendix B: Slab with a local constituent permittivity: illustrative phenomenological model 


In order to illustrate the general nature of the negative static permittivity and, consequently, of the violation of 
KK relations in Q2D systems, in this section we consider a simple model of a slab of the thickness a comprised of a 
uniform medium with the local permittivity e(w), as shown in Fig. 10. Let the externally applied potential be 

:(>ext(r,a;) = e*('i-—(Bl) 

where q lies in the xy plane. Then the total potential in the regions I, II, and III can be written as 


(j){z) = 


I l + Ale-«^ 

II 1 -b A2e-«^ -b B2e^^ 

III l-bB3e«^ 


(B2) 














7 


I 6=1 



z 



FIG. 10. A slab with the local permittivity £(uj) and the thickness a surrounded by vacuum. 


where Ai, A 2 , B 2 , and B 3 are constants to be found from the boundary conditions 

£ 

1 + = -+ + B2e“«“/^ 

£ 


Bse 


— qaf2 _ 


= £ 




where the first two are due to the continuity of 4 >{z) at the interfaces and the last two to the continuity 
z-component of the displacement vector = —£(z) . The solution of Eqs. (B3) is 


A2 — B2 — 


e — 1 




^ £ 1 + £“9 + £(e“'? - 1) ’ 

Ai = S3 =£(l-e“«)A 2 . 

The induced charge-density can be found as 

- 9') ['^(^) - '^ext(^)] . 

Therefore, using Eqs. (B2) and 0ext(-z) = 1, we have 

® (i ■ ^) ® (i+^)+■ i) [(^1 ■ 

+d (z -f [(^3 - B 2 ) e-«“/2 + } , 

where ©(x) is the Heaviside step-function. Then 

J x(z, z')dzdz' = J p{z)dz = ~ ^ Aie~'^°'^^ + 2 A 2 sinh {qa/ 2 ) 

where Eqs. (B4) were used. Einally, using Eq. (1), we have 

1 




^=1+^ 

e{q,u;) ^ 2 


£(w) 


- 1 


Ai(q,a;)e -|- 2^2(9, w) sinh (ga/2), 


where we have restored the explicit arguments of the functions. 

In Fig. 11 we plot the permittivity of Eq. (B8) for a metallic slab with the Drude constituent permittivity 


OJ" 


(w -I- iqY 


(B3) 


of the 

(B4) 


(B5) 


(B6) 


(B7) 


(B8) 


£(w) = 1 


(B9) 



















aq = 0.2 



FIG. 11. The permittivity of the metallic slab model below (left) and above (right) the critical wave-number. 


with rj/ujp = 0.067. Similar to the first-principles calculations for Q2D crystals, the model system yields the negative 
static permittivity at lager wave-vectors (right panel of Fig. 11). Finally, by the direct substitution of Eqs. (B4) and 
(B9) into Eq. (B8), it is easy to show that, in the limit of the 2D electron gas (a —)• 0, oj^a = innsDCi —>■ 47rn2D), the 
permittivity (B8) reduces to 


e{q,u}) 1 


27rrt2Dg 
(w -I- ir])'^ ’ 


(BIO) 


which is a standard result for the 2D electron gas in the long-wave limit, with which, the negative static permittivity 
is not, of course, possible. 


Appendix C: Beyond RPA: Adiabatic time-dependent local-density approximation 



FIG. 12. Gomparison of the permittivities of a single-layer graphene within RPA (black and red lines) and ATDLDA (blue and 
green lines). Left: q < qc. Right: q > qc- 


In order to go beyond RPA, we have conducted calculations using the xc kernel fxc at the level of ATDLDA. Results 
are presented in Fig. 12, showing no significant difference compared with RPA. It must be noted that for TDDFT as 
applied to Q2D crystals, ATDLDA is the current state-of-the-art in accounting for the dynamic xc effects. Indeed, 
more elaborate kernels [27, 28], developed for 3D crystals, are not applicable to the Q2D case. 
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Appendix D: Distribution of the charge-density and the total potential in the ^-direction 
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FIG. 13. The 2 -distribution of the charge-density (upper panel) and total potential (lower panel) at two values of the wave- 
vector. The static (ui = 0), uniform in the 2 -direction, unity-amplitude external potential is applied. The averaging in the 
ly-plane has been performed. 


In Fig. 13 we plot the charge-density p{q,z) induced in graphene and the corresponding total potential 4>{q,z) in 
response to the static, uniform in the z-direction, unity-amplitude external potential for the two values of the in-plane 
wave-vector, below and above the critical value qc- An important feature seen from the lower panel is that, although a 
strong screening occurs inside the graphene layer, the direction of the total field does not change sign even for q > q^. 
This prevents us from directly relating the negative static permittivity in Q2D crystals to the 2D superconductivity 
in these systems. 










